실습: 경사 충격파 관계식#

강좌: 수치해석 프로젝트

이론#

수직충격파#

초음속 유동이 물체이 부딛히거나 유동의 방향이 꺽이는 경우 충격파 또는 팽창파가 발생한다. 수직 방향으로 충격파가 발생하는 경우 충격파 전/후 마하수는 다음과 같은 Prandtl 관계를 만족한다.

\[ M_1^{*} M_2^{*} = 1 \]

여기서 \(M^*\) 는 특성 마하수로 Sonic point일 때 음속으로 마하수 \(M\) 과는 아래 관계를 만족한다.

\[ M^2 = \frac{2}{(\gamma+1)/M^{*2} - (\gamma-1)}, M^{*2} = \frac{(\gamma+1) M^2}{2 + (\gamma-1)M^2} \]

여기서 \(\gamma\) 는 비열비로 일반 공기는 1.4이다.

이 식을 이용하면 Prandtl 관계식은 다음과 같이 정리할 수 있다.

\[ M^2_2 = \frac{1 + [(\gamma-1)/2]M_1^2}{\gamma M_1^2 - (\gamma-1)/2} \]

충격파 전/후 물성치 변화는 다음과 같다.

\[ \frac{\rho_2}{\rho_1} = \frac{(\gamma+1)M_1^2}{2 + (\gamma-1)M_1^2} \]
\[ \frac{p_2}{p_1} = 1 + \frac{2\gamma}{\gamma+1} (M_1^2-1) \]

충격파 전/후에는 전압력 손실이 발생한다. 충격파를 제외한 영역은 등엔트로피 유동으로 가정할 수 있다. 즉 마하수가 \(M\) 인 유동의 정압 (Static Pressure, \(p\)) 대비 전압력 (Total Pressure, \(p_0\)) 는 다음 관계를 갖는다.

\[ \frac{p_0}{p} = \left ( 1 + \frac{\gamma-1}{2} M^2 \right )^{\gamma/(\gamma-1)}. \]

전압력 손실은 충격파 전/후의 전압력을 각각 구한 후 그 감소 비율이다.

경사충격파#

쐐기 각도 \(\theta\) 인 경우 충격파의 파각은 \(\beta\)로 발생한다. 이 때 충격파에 수직인 성분은 수직충격파 관계식을 따르고, 수평 방향 속도는 같다. 즉

\[ Vn_1 = V_1 \sin(\beta), Vn_2 = V_2 \sin(\beta - \theta) \]
https://www.researchgate.net/profile/Buraq-Al-Mosawi/publication/332299135/figure/fig1/AS:745770600902656@1554816975212/Fig-1-Geometry-of-flow-through-an-oblique-shock-wave.ppm

Fig. 3 Shock Wave#

\(\theta-\beta-M\) 에 대해서는 아래 관계식을 만족한다.

\[ \tan \theta = 2 \cot \beta \frac{M_1^2 \sin^2\beta - 1}{M_1^2 (\gamma+\cos2\beta)+2} \]
%matplotlib inline
from matplotlib import pyplot as plt

import numpy as np

plt.style.use('ggplot')
plt.rcParams['figure.dpi'] = 150
def tbm(beta, M, gamma=1.4):
    """
    Theta-Beta-M relation
    
    Parameters
    ----------
    beta : float
        Angle of wave (rad)
    M : float
        Mach number
    gamma : float
        specific heat ratio
        
    Returns
    -------
    theta : float
        Angle of wedge (rad)
    """
    # tangent(theta)
    tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
    
    # tan^{-1}(theta)
    return np.arctan(tangent)
M = 1.1
beta = np.deg2rad(10)
theta = np.rad2deg(tbm(beta, M))
print('Theta={} @ M={}, beta={}'.format(theta, M, beta))

beta = np.deg2rad(70)
theta = np.rad2deg(tbm(beta, M))
print('Theta={} @ M={}, beta={}'.format(theta, M, beta))
Theta=-66.15222865219691 @ M=1.1, beta=0.17453292519943295
Theta=1.0317298907919357 @ M=1.1, beta=1.2217304763960306
from scipy import optimize


def tbm_curve(M, gamma=1.4, n=31):
    # Construct equation for beta
    f = lambda x : tbm(x, M, gamma)
    
    # Compute root of tbm w.r.t beta
    sol = optimize.root_scalar(f, bracket=[0, np.pi/2])
    beta_min = sol.root

    # Divide betas on (beta_min, 90)
    beta = np.linspace(beta_min, np.pi/2, n)
    
    # Compute thetas
    theta = tbm(beta, M)

    return theta, beta
# Plot for Mach numbers
for M in [1.1, 1.2, 1.4, 1.6, 1.8, 2, 2.2, 2.5, 3, 4, 6, 15]:
    # Get theta-beta curve (in degree)
    theta, beta = np.rad2deg(tbm_curve(M))
    
    # Plot each curve
    plt.plot(theta, beta, color='gray')
    
    # Get maximum theta and beta 
    idx = theta.argmax()
    theta_max = theta[idx]
    beta_loc = beta[idx]
    
    # Add text
    plt.text(theta_max, beta_loc, "M={}".format(M), fontsize='x-small')
    
# Add labels using Latex symbols
plt.xlabel(r"$\theta$")
plt.ylabel(r"$\beta$")
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
  tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
  tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
  tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
  tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
  tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
  tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
  tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
  tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
  tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
  tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
  tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
  tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
Text(0, 0.5, '$\\beta$')
../_images/e3033718ccd13f02c0528963706bfa426355dd953833fa4c397a1a69a2d15ebc.png

실습#

수직충격파#

  • 마하수에 따른 전압력을 계산하는 코드를 작성하시오.

  • 충격파 전 마하수 M1 일 때 충격파 후 마하수 및 밀도 및 압력비를 구하는 코드를 작성하시오

def p0_p(M, gamma=1.4):
    #...
    return ...


def normal_shock(M1, gamma=1.4):
    # ...
    return M2, rho2, p2, p0_loss

경사충격파#

  • 충격파 전 마하수 M1, 쐐기 각도가 theta 일 때 \(\theta-\beta-M\) 비선형 방정식을 해석하시오.

    • \(f(\beta, M) \rightarrow \tan(\theta\)) 를 구하는 코드를 만드시오.

    • \(\beta\)가 달라질 때 \(-\tan(\theta)\) 를 최소화하는 조건과 그 때 \(\theta\)를 구하시오.

    • \(\beta\) 를 미지수 x 라 하고, 주어진 \(M, \theta\) 에 대해서 \(f(x, M) - \tan(\theta)=0\) 인 방정식의 해를 수치적으로 구하시오.

  • 충격파 후 마하수, 밀도, 압력비 및 파각을 구하는 코드를 작성하시오.

    • Hint) 파각에 수직방향 마하수 \(M_{n1}, M_{n2}\)

    \[ M_{n1} = \frac{M_1}{\sin(\beta)}, M_{n2} = \frac{M_2}{\sin(\beta-\theta)} \]
def tangent_theta(beta, M, gamma=1.4):
    # Calculuate theta-beta-M relation
    return t_theta


def theta_max(M, gamma=1.4):
    # function to compute -tangent_theta(x, M, gamma=1.4) for given M
    
    # minimize the function with very small x0
    
    return theta


def beta_week(M, theta, gamma=1.4):
    # make equation f(x,M) - tan(theta) =0
    
    # root finding
    
    return beta


def oblique_shock(M1, theta, gamma=1.4):
    # Solve theta-beta-M
    
    # Compute Mn1
    
    # Solve normal shock for Mn1
    
    # Calculate M2
    
    return M2, rho2, p2, p0_loss beta